#################################
# ANALYSIS: EFFECTS OF INDIRECT RULE
# Effect of v33 x colonizer on current TPIs
# 
# Part of replication of:
# Traditional Institutions in Africa, Past and Present
#
# Political Science Research and Methods
#
# By Clara Neupert-Wentz and Carl Müller-Crepon, 2023
#
##########################################



# PREPARE ANALYSIS ##############

## Define main data
data <- readRDS("data/main_data.rds")

## Main outcome
main.depvar <- c(`TPI Index` = "tpi_pc1")

all.depvar <- plot.vars <- c( "tpi_pc1", "tni_level_ord", 
                              "leader_max",
                              "inst_mean","leader_mean", "func_mean", "ties_mean")
all.depvar.labs <- plot.labs <- c("TPI Index", "TPI Level",
                                   "Max Leader", 
                                  "Institution Index","Leader Index",  "Functions Index", "State-ties Index")
names(all.depvar) <- all.depvar.labs


## Sort version -- 3 indices
short.depvar  <- c( "tpi_pc1",  "tpi_centr_idx", "tpi_funct_idx")
short.depvar.labs <- c("TPI Index", "Political Centralization", "Functional Differentiation")
names(short.depvar) <- short.depvar.labs


## Main Treatment
main.treat <- c(`Jurisdictional Hierarchy (v33)` = "v33.num")


## Define main covariates

### Baseline
base.controls <- c("log(1 + pop.1880)", "log(poly.area)",
                   "log(coast.dist)", "log(river.dist)")
base.controls.lab <- c("Population","Area", "Distance to coast", "Distance to nav. river")
names(base.controls) <- base.controls.lab

### Ethnic vars from Murdock
ethnic.controls <- paste0("v",c(4,5,28),".num")
ethnic.controls.lab <- c("Reliance on agriculture",
                         "Reliance on pastoralism", "Intensity of agriculture")
names(ethnic.controls) <- ethnic.controls.lab

### Nature
nature.controls <- c("medianaltitude","medianslope","temperaturemean","evapotranspiration",  
                     "precipitation", "ppetratio","suit", "max.cash.crop",
                     "malaria.tsuit.max", "tsetse.max")
nature.controls.lab <- c("Altitude","Ruggedness","Temperature","Evapotranspiration",  
                         "Precipitation", "Precipitation/Evapotr.","Agr. suitability", "Cash crop suitability",
                         "Malaria environment", "Tsetse environment")
names(nature.controls) <- nature.controls.lab

### Combine into main three specifications
main.cov.ls <- list(c(),
                    c(base.controls),
                    c(base.controls,ethnic.controls),
                    c(base.controls,ethnic.controls,nature.controls))

## Cluster of SEs?
cluster.var <- "gid" ## Murdock group (across borders)


## Main fixed effects
main.fe <- "iso3c" ## Country

# Interaction specification

## Treatments
treat.int <- c("british", paste0(main.treat, c("_british", "_french")))

## Covariates (full interactions)
int.cov.ls <- lapply(main.cov.ls, function(cv){if(length(cv) > 0){c(cv, paste0("I(british*", cv,")"))} else {cv}})

## Stargazer Output

### What to keep?
star.keep <- 1:3
names(star.keep) <- c("British",
                      "Juris. hier. (v33) $\\times$ British", 
                      "Juris. hier. (v33) $\\times$ French" )


### Output type
star.out <- "latex" ### At some point, this will become nicer

### Notes
if(star.out == "latex"){
  latex.notes.add <- function(width = .95, cluster = "ethnic group level"){
    paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS models. Standard errors are clustered on the ethnic group level. Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}
  latex.notes.add.int <- function(width = .95, cluster = "ethnic group level"){
    paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS models. Standard errors are clustered on the ethnic group level.
           Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}
  
} else {
  latex.notes.add <- function(){""}
}




# DESCRIPTIVE PLOT ###
# Figure 2, main text


# Data
plot.df <- data
plot.df$colonizer <- "Full Sample"
plot.df <- rbind(plot.df,
                 data[data$colonizer %in% c("British", "French"),])
plot.df$colonizer <- factor(plot.df$colonizer, 
                            levels = c("Full Sample", "British", "French"),
                            labels = c("Full Sample", "Former British Rule", "Former French Rule"),
                            ordered = T)

# Summary
sum.df <- summarySE(data=plot.df, measurevar = "tpi_pc1", groupvars= c("colonizer", "v33.num"), na.rm=T,
                    conf.interval=.95, .drop=TRUE)

# Plot Scatterplot matrix
g <- ggplot(plot.df, aes(x = v33.num, y = tpi_pc1)) +
  geom_jitter(width = .2, height = .05, alpha = .1) + 
  geom_smooth(method = "lm", color = "blue") +
  geom_errorbar(data = sum.df, 
                aes(ymin=tpi_pc1-ci, ymax=tpi_pc1+ci), width=.15, lwd = 1.3,
                color = "white") +
  geom_point(data = sum.df, size = 1.5, col = "white") +
  geom_point(data = sum.df, size = .75) +
  geom_errorbar(data = sum.df, 
                aes(ymin=tpi_pc1-ci, ymax=tpi_pc1+ci), width=.1, lwd = .7) +
  ylab("TPI Index (0-1)") +
  xlab("Jurisdictional hierarchy (Murdock's v33)") +
  facet_wrap(~ colonizer, nrow = 1) +
  theme_minimal()  + theme(panel.spacing = unit(1, "lines")) +
  NULL

png(file.path(fig.path, paste0("scatter_indrule.png")), width = 6, height = 2.5, unit = "in", res = 400)
print(g)
dev.off()



# SHORT OUTCOMES #################
# Figure 5, main text

# Estimate
model.ls <-  unlist(lapply(short.depvar, function(dv){
  spec.ls <- lapply(main.cov.ls, function(cv){
    make_form(dv = dv,
              expl = c(treat.int, cv),
              fe = main.fe, iv = "0", se = cluster.var)
  })
  
  lapply(spec.ls, function(spec){
    felm(spec, data = data)
  })
}), recursive = F)


# Get coefficients
coef.ls <- extract_coef(model.ls)

# Prepare plot
plot.df <- do.call(rbind, lapply(coef.ls, function(c){
  data.frame(outcome = c$dv,
             treat = c("British", "French", "Difference"), 
             coef = c(c$coef[treat.int[2], 1], c$coef[treat.int[3], 1], 
                      c$coef[treat.int[3], 1] - c$coef[treat.int[2], 1]),
             se = c(c$clustervcv[treat.int[2], treat.int[2]],
                    c$clustervcv[treat.int[3], treat.int[3]], 
                    c$clustervcv[treat.int[2], 
                                 treat.int[2]] + 
                      c$clustervcv[treat.int[3], treat.int[3]] - 
                      2*c$clustervcv[treat.int[3], treat.int[2]])^.5)
}))
plot.df$spec <- rep(rep(1:length(int.cov.ls), length(short.depvar)), each = 3)
plot.df$treat <- factor(plot.df$treat, levels = unique(plot.df$treat), ordered = T)
plot.df$treat.num <- as.numeric(plot.df$treat)
plot.df$outcome_label <- factor(rep(short.depvar.labs, each = length(int.cov.ls) * 3), 
                                levels = short.depvar.labs, ordered = T)

# Plot


## Prepare
g <- ggplot(plot.df, aes(y = coef, x = treat.num + ((spec-2.5)/6), group = treat, color = treat)) + 
  geom_point(size = .75) + 
  geom_errorbar(aes(ymin = coef - 1.96*se, 
                    ymax = coef + 1.96*se), 
                width=0, lwd = .5) +
  scale_color_manual(values=c("red", "blue", "black")) +
  geom_hline(yintercept = 0, color = "darkgrey", lty = 2) +  
  scale_x_continuous(breaks=c(1:3), labels = c("British", "French", "Diff."),
                     limits = c(0,4))  + 
  theme_minimal() +
  xlab("Treatment x Specification") + ylab("Marginal effect of v33 on outcome") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.35, hjust = 1),
        legend.position = "none") +
  facet_wrap(~ outcome_label, nrow = 1) + 
  labs(title = "Long-term effect of British vs. French rule",
       subtitle = "Each cluster of coefficients plots the results from the four main specifications.") +
  NULL

## Print
png(file.path(fig.path, paste0("indrule_shortoutcomes.png")), width = 6.5, height = 3, unit = "in", res = 400)
par(mar = c(0,0,0,0))
print(g)
dev.off()



# ALL OUTCOMES #################
# Figure 6, main text

# Estimate
model.ls <-  unlist(lapply(all.depvar, function(dv){
  spec.ls <- lapply(main.cov.ls, function(cv){
    make_form(dv = dv,
              expl = c(treat.int, cv),
              fe = main.fe, iv = "0", se = cluster.var)
  })
  
  lapply(spec.ls, function(spec){
    felm(spec, data = data)
  })
}), recursive = F)


# Get coefficients
coef.ls <- extract_coef(model.ls)

# Prepare plot
plot.df <- do.call(rbind, lapply(coef.ls, function(c){
  data.frame(outcome = c$dv,
             treat = c("British", "French", "Difference"), 
             coef = c(c$coef[treat.int[2], 1], c$coef[treat.int[3], 1], c$coef[treat.int[3], 1] - c$coef[treat.int[2], 1]),
             se = c(c$clustervcv[treat.int[2], treat.int[2]],
                    c$clustervcv[treat.int[3], treat.int[3]], 
                    c$clustervcv[treat.int[2], treat.int[2]] + c$clustervcv[treat.int[3], treat.int[3]] - 2*c$clustervcv[treat.int[3], treat.int[2]])^.5)
}))
plot.df$spec <- rep(rep(1:length(int.cov.ls), length(plot.vars)), each = 3)
plot.df$treat <- factor(plot.df$treat, levels = unique(plot.df$treat), ordered = T)
plot.df$treat.num <- as.numeric(plot.df$treat)
plot.df$outcome_label <- factor(rep(all.depvar.labs, each = length(int.cov.ls) * 3), levels = all.depvar.labs, ordered = T)

# Plot


## Prepare
g <- ggplot(plot.df, aes(y = coef, x = treat.num + ((spec-2.5)/6), group = treat, color = treat)) + 
  geom_point(size = .75) + 
  geom_errorbar(aes(ymin = coef - 1.96*se, 
                    ymax = coef + 1.96*se), 
                width=0, lwd = .5) +
  scale_color_manual(values=c("red", "blue", "black")) +
  geom_hline(yintercept = 0, color = "darkgrey", lty = 2) +  
  scale_x_continuous(breaks=c(1:3), labels = c("British", "French", "Diff."),
                     limits = c(0,4))  + 
  theme_minimal() +
  xlab("Treatment x Specification") + ylab("Marginal effect of v33 on outcome") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.35, hjust = 1),
        legend.position = "none") +
  facet_wrap(~ outcome_label, nrow = 1) + 
  labs(title = "Long-term effect of British vs. French rule",
       subtitle = "Each cluster of coefficients plots the results from the four main specifications.") +
  NULL

## Print
png(file.path(fig.path, paste0("indrule_alloutcomes.png")), width = 7, height = 3, unit = "in", res = 400)
par(mar = c(0,0,0,0))
print(g)
dev.off()

########################################
# APPENDIX FIGURES #####################

# DESCRIPTIVE PLOT 3 OUTCOMES ###
# Appendix Figure A2

# Data
plot.df <- do.call(rbind, lapply(seq_along(short.depvar), function(i){
  plot.df <- data
  plot.df$colonizer <- "Full Sample"
  plot.df <- rbind(plot.df,
                   data[data$colonizer %in% c("British", "French"),])
  plot.df$colonizer <- factor(plot.df$colonizer, 
                              levels = c("Full Sample", "British", "French"),
                              labels = c("Full Sample", "Former British Rule", "Former French Rule"),
                              ordered = T)
  plot.df$xvar <- plot.df[, short.depvar[i]]
  plot.df$xlab <- short.depvar.labs[i]
  plot.df
}))
plot.df$xlab <- factor(plot.df$xlab, levels = unique(plot.df$xlab), ordered = T)

# Summary
sum.df <- summarySE(data=plot.df, measurevar = "xvar", 
                    groupvars= c("colonizer", "xlab", "v33.num"), na.rm=T,
                    conf.interval=.95, .drop=TRUE)

# Plot Scatterplot matrix
g <- ggplot(plot.df, aes(x = v33.num, y = xvar)) +
  geom_jitter(width = .2, height = .05, alpha = .1) + 
  geom_smooth(method = "lm", color = "blue") +
  geom_errorbar(data = sum.df, 
                aes(ymin=xvar-ci, ymax=xvar+ci), width=.15, lwd = 1.3,
                color = "white") +
  geom_point(data = sum.df, size = 1.5, col = "white") +
  geom_point(data = sum.df, size = .75) +
  geom_errorbar(data = sum.df, 
                aes(ymin=xvar-ci, ymax=xvar+ci), width=.1, lwd = .7) +
  ylab(NULL) +
  xlab("Jurisdictional hierarchy (Murdock's v33)") +
  facet_grid(xlab ~ colonizer, switch = "y") +
  theme_minimal()  + 
  theme(panel.spacing = unit(1, "lines"),
        strip.placement = "outside",
        strip.text = element_text(size = 11)) +
  NULL

png(file.path(fig.path, paste0("scatter_indrule_3out.png")), 
    width = 6, height = 7, unit = "in", res = 400)
print(g)
dev.off()




#############################
# APPENDIX TABLES ###########

# BASELINE MODEL ############
# Appendix Table A4

# Stub
stub <- "indrule_main"

# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = main.depvar,
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Jurisdictional hierarchy and current TPIs in former British and French colonies",
            dep.var.caption = names(main.depvar),dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)




# SPLIT TPI INDEX ############
# Appendix Table A7

# Stub
stub <- "indrule_splittpiidx"


# Setup equations -- Fixed effects
spec.ls <- unlist(lapply(c("tpi_centr_idx","tpi_funct_idx"), function(dv){
  lapply(main.cov.ls[c(1,4)], function(cv){
    make_form(dv = dv,
              expl = c(treat.int, cv),
              fe = main.fe, iv = "0", se = cluster.var)
  })
}), recursive = F)

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", rep(c("no", "yes"), 2)),
                  latex.addline("Ethnic covariates", rep(c("no",  "yes"), 2)),
                  latex.addline("Nature covariates", rep(c("no", "yes"), 2)),
                  latex.addline("Country (2016) FEs", rep(c("yes", "yes"), 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Persistence in former British and French colonies: Splitting the TPI Index",
            column.labels = c("Centralization Idx", "Functions Idx"),
            column.separate = c(2,2),
            dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)







# NO FIXED EFFECTS ################
# Appendix Table A9

# Stub
stub <- "indrule_nofe"

# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = main.depvar,
            expl = c(treat.int, cv),
            fe = "0", iv = "0", se = cluster.var)
})

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(
  latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){
    paste0(round(x$Estimate, 3), 
           ifelse(x$pvalue < .01, "^{***}",
                  ifelse(x$pvalue < .05, "^{**}",
                         ifelse(x$pvalue < .1, "^{*}",
                                "" ) ) ))
  }))),
  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
  latex.addline("Country (2016) FEs", c("no", "no", "no", "no")),
  latex.mean.dv(model.ls)
)
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Jurisdictional hierarchy and current TPIs in former British and French colonies",
            dep.var.caption = names(main.depvar),dep.var.labels.include = FALSE,
            keep= 1:4,
            covariate.labels = c("Constant", names(star.keep)),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)



# FULL INTERACTIVE MODEL ############
# Appendix Table A10
# Stub
stub <- "indrule_fullint"

# Setup equations -- Fixed effects
spec.ls <- lapply(int.cov.ls, function(cv){
  make_form(dv = main.depvar,
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline $\\times$ British covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic $\\times$ British covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature $\\times$ British covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Former British and French colonies: Full interactions",
            dep.var.caption = names(main.depvar),dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)



# NO STATE TIES ############
# Appendix Table A11

# Stub
stub <- "indrule_nostateties"


# Setup equations -- Fixed effects
spec.ls <- unlist(lapply(paste0(short.depvar, "_nst")[c(1,3)], function(dv){
  lapply(main.cov.ls[c(1,4)], function(cv){
  make_form(dv = dv,
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})}), recursive = FALSE)

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences

diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no",  "yes", "no",  "yes")),
                  latex.addline("Ethnic covariates", c("no", "yes", "no", "yes")),
                  latex.addline("Nature covariates", c("no",  "yes","no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes","yes", "yes","yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="TPI and Function Index w/out state ties",
            dep.var.labels.include = FALSE,
            column.labels = short.depvar.labs[c(1,3)], 
            column.separate = c(2,2),
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)






# IMPORTANCE ############
# Appendix Table A12

# Stub
stub <- "indrule_importance"


# Setup equations -- Fixed effects
spec.ls <- unlist(lapply(c("I(imp_day/5)","I(imp_nat/5)"), function(dv){
  lapply(main.cov.ls[c(1,4)], function(cv){
    make_form(dv = dv,
              expl = c(treat.int, cv),
              fe = main.fe, iv = "0", se = cluster.var)
  })
}), recursive = F)

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", rep(c("no", "yes"), 2)),
                  latex.addline("Ethnic covariates", rep(c("no",  "yes"), 2)),
                  latex.addline("Nature covariates", rep(c("no", "yes"), 2)),
                  latex.addline("Country (2016) FEs", rep(c("yes", "yes"), 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Importance of traditional institutions",
            column.labels = c("Daily life", "National politics"),
            column.separate = c(2,2),
            dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)











# MATCHING  TPI - MURDOCK MODEL ############
# Appendix Table A13

# Stub
stub <- "indrule_matched"

# Matched dummy
data$matched <- ifelse(data$group == "NA", 0, 1)

# Setup equations -- Fixed effects
spec.ls <- c(list( make_form(dv = "matched",
                             expl = c(treat.int),
                             fe = "0", iv = "0", se = cluster.var)),
             lapply(main.cov.ls, function(cv){
  make_form(dv = "matched",
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
}))

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("no", "yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Selection stage 1: Selection into sample of Murdock groups matched to TradGovGroups",
            dep.var.caption = "Match of Murdock group with TradGov group",
            dep.var.labels.include = FALSE,
            keep= 1:4, covariate.labels = c("Constant", names(star.keep)),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)



# ANY TPI MODEL ############
# Appendix Table A14

# Stub
stub <- "indrule_anytpi"


# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = "any_tpi",
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Selection stage 2: Selection into non-missing data on traditional institutions",
            dep.var.caption = "Any TPI",dep.var.labels.include = FALSE,
            keep= star.keep, covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)


# FULL SELECTION MODEL ############
# Appendix Table A15

# Stub
stub <- "indrule_select"


## Selection dummy
data$select <- as.numeric(data$matched == 1 & data$any_tpi > 0)


# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = "select",
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Selection stage 1 + 2: Selection of Murdock groups into final sample",
            dep.var.caption = "Murdock group in final sample",
            dep.var.labels.include = FALSE,
            keep= star.keep, covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)


# REWEIGHTING TO ACCOUNT FOR SELECTION ##########
# Appendix Table A16

# Stub
stub <- "indrule_selectreweight"

# selection

## Selection model for French subset
sel.m <- glm(as.formula(paste("select ~ ", 
                              paste(c("v33.num",main.cov.ls[[length(main.cov.ls)]]), collapse = "+"))),
             data = data[data$french == 1,], family  = "binomial")
summary(sel.m)

## Selection probability
sel.prob <- predict.glm(sel.m, newdata = data, type = "response")

# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = main.depvar,
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate with weights from above
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data[!is.na(sel.prob),], 
       weights = 1 / sel.prob[!is.na(sel.prob)])
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Former British vs. French colonies: Reweighting by estimated `selection' probability",
            dep.var.caption = names(main.depvar),dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)



# RECODE MISSINGS TO ZERO ##########
# Appendix Table A17

# Stub
stub <- "indrule_rectozero"

# Recode
data$tpi_pc1_rec <- ifelse(data$select == 1, data$tpi_pc1, 0)


# Setup equations -- Fixed effects
spec.ls <- lapply(main.cov.ls, function(cv){
  make_form(dv = "tpi_pc1_rec",
            expl = c(treat.int, cv),
            fe = main.fe, iv = "0", se = cluster.var)
})

# Estimate with weights from above
model.ls <-  lapply(spec.ls, function(spec){
  felm(spec, data = data)
})

# Calculate differences
diff.ls <- lapply(model.ls, function(m){
  vars <- treat.int[2:3]
  dm <- deltaMethod(m, paste(vars, collapse = " - "))
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(vars, collapse = " - "))), type = "robust")
  dm$pvalue <- wt["p"]
  dm
})

# Prepare table
add.lines <- list(latex.addline_nomc("British-French Diff.:", unlist(lapply(diff.ls, function(x){paste0(round(x$Estimate, 3), 
                                                                                                        ifelse(x$pvalue < .01, "^{***}",
                                                                                                               ifelse(x$pvalue < .05, "^{**}",
                                                                                                                      ifelse(x$pvalue < .1, "^{*}",
                                                                                                                             "" ) ) ))}))),
                  latex.addline_nomc("", unlist(lapply(diff.ls, function(x) paste0("(",round(x$SE, 3), ")")))),
                  latex.addline("Baseline covariates", c("no", "yes", "yes", "yes")),
                  latex.addline("Ethnic covariates", c("no", "no", "yes", "yes")),
                  latex.addline("Nature covariates", c("no", "no", "no", "yes")),
                  latex.addline("Country (2016) FEs", c("yes", "yes", "yes", "yes")),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# Print
fileConn<-file(file.path(tab.path, paste0(stub, ".tex")))
writeLines(
  stargazer(model.ls,
            title="Former British vs. French colonies: Recoding missing groups to zero",
            dep.var.caption = names(main.depvar),dep.var.labels.include = FALSE,
            keep= star.keep,  covariate.labels = names(star.keep),
            notes.align = "l",label=paste0(stub),align =T,
            add.lines = add.lines, 
            digits = 3, intercept.top = T,intercept.bottom = F,
            omit.stat = c("rsq","res.dev","ser"),column.sep.width = "8pt",
            notes = latex.notes.add(), notes.label = "", notes.append = F,
            type  = star.out), 
  fileConn)
close(fileConn)






